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The pion form factor is computed using quenched twisted mass QCD and the 
GMRES-DR matrix inverter. The momentum averaging procedure of Frezzotti and 
Rossi is used to remove leading lattice spacing artifacts, and numerical results for 
the form factor show the expected improvement with respect to the standard Wilson 
action. Although some matrix inverters are known to fail when applied to twisted 
mass QCD, GMRES-DR is found to be a viable and powerful option. Results ob- 
tained for the pion form factor are consistent with the published results from other 
0{a) improved actions and are also consistent with the available experimental data. 



I. OVERVIEW AND DISCUSSION 



The Wilson lattice action for fermionsjlj does not have exact chiral symmetry, and 
thus it provides no lower bound for the norm of Dirac matrix eigenvalues. As a practi- 
cal consequence, numerical simulations encounter "exceptional configurations" (for which 
the Dirac matrix cannot be inverted) with increasing regularity as the quark mass is re- 
duced. Symanzik improvement of the Wilson action by adding the Sheikholeslami-Wohlert 
term with a properly tuned coefficient [3] , exacerbates the problem. 

Fortunately, the addition of a twisted mass term provides a lower bound for the eigen- 
values thereby eliminating exceptional configurations. Furthermore, tuning of the 
Sheikholeslami-Wohlert coefficient is no longer needed for the removal of linear lattice spacing 
errors. Instead, many quantities are automatically improved simply by adding the twisted 
mass term to the basic Wilson action and setting the hopping parameter to its critical 
value (i.e. k = in Eq. (fT^ below), while other quantities become improved by averaging 
over equal and opposite momenta. 0, 0] The present work provides an explicit numerical 
verification of the averaging technique for improvement in twisted mass QCD (tmQCD). 

Our discussion will focus on the pion electromagnetic form factor. Due to the simple 
valence quark structure of the pion and a firm theoretical knowledge of its form factor in 
both the and — > 00 limits, the pion form factor is a preferred place to study 
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the transition between perturbative and nonperturbative QCD. Experimental studies are 
ongoing at Jefferson Lab, and tlieoretical modelling is also continuing. ^ 



Initial studies of the pion form factor using lattice QCD occurred some time ago,[2l,IS| and 
new lattice initiatives have arisen recently for Wilson, Sheikoleslami-Wohlert, and domain 
wall actions. 0, 1^ JJ, | Preliminary results from our tmQCD study were presented in Ref. . 
In contrast to other form factors, the pion form factor receives no contributions from the 
vector current attaching to a nonvalence quark (so-called "disconnected diagrams"), and this 
feature reduces the lattice QCD cost considerably. There could still be contributions from 
sea quarks that do not interact directly with the external vector current, and these have 
been considered in Ref. where dynamical configurations were to used obtain the pion 
form factor. All other studies to date, including the present one, have used the quenched 
approximation and thereby omitted all nonvalence quarks. 

In the remainder of this article, we report on our use of quenched tmQCD to compute the 
pion form factor. Two quark masses corresponding to pion masses near 470 MeV and 660 
MeV, as well as a variety of momentum transfer values satisfying GeV^ < ^ 5 GeV^ 
have been considered, all at a lattice spacing of 0.10 fm. A comparison to existing lattice 
results clearly shows the improvement expected for tmQCD, since the momentum-averaged 
tmQCD data agree with results from other improved actions and differ from unimproved 
Wilson results at this same lattice spacing. Interestingly, even before momentum averaging 
the tmQCD data are closer to improved action results than to unimproved Wilson results 
for this particular observable, despite the fact that the pion form factor technically requires 
momentum averaging to exactly remove the linear lattice spacing errors. 

To determine the renormalization factors that appear in the pion form factor correlation 
function, and to compare with the predictions of vector meson dominance, we also study two- 
point pseudoscalar and vector correlators with nonzero momenta. The associated dispersion 
relations are compared to continuum expectations as another means of exploring lattice 
spacing artifacts. 

One of the technical issues that arises in tmQCD simulations is the failure of some 
standard matrix inversion algorithms. Alternative algorithms are being used and evaluated 
by various authors. The present work makes use of the GMRES-DR algorithm [isjl and 
concludes that it performs well for tmQCD. Some details are presented in Section HVl 

This initial exploration of the pion form factor with tmQCD leads to optimism that 
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future lattice tmQCD studies, perhaps with smeared operators and increased statistics on 
larger lattices, can reach smaller quark masses with greater precision. Our present results 
are consistent with vector meson dominance and with experiment. More generally, the 
present work underscores the value of lattice tmQCD itself as a practical tool for hadron 
phenomenology. 

II. CORRELATION FUNCTIONS 

The electromagnetic form factor of a charged pion is defined by 

(vr+(p7)b.(0)k+(pl)) = F{Q'){p,+pf), (1) 

where j^{0) is a conserved vector current evaluated at the spacetime origin, pi and pf are the 
initial and final pion (Euclidean) 4-momenta respectively, Pi and pf are the corresponding 
3-momenta, and = {pj — Pi)^ is the 4-momentum transfer. To compute this matrix 
element on a spacetime lattice, one can use the three-point correlator displayed in Fig. ^ 
A source with pion quantum numbers is placed at Xj, a sink at a;/, and a vector current is 
inserted at x. Given an interpolating field operator with the quantum numbers of the 
charged pion, one can extract the form factor from the following three-point correlator: 

r.^.(t.,t,t/,pl,P/) = E (T^^^-^-^^er^^-^^-^^ (O|0(x;)j^(x)0t(a;,)|o) . (2) 

Note that we are working in units of lattice spacing throughout this discussion. In this work 
we choose to be the local operator 

<\){x) = d{x)'j5u{x) (3) 

where u{x) and d{x) are the up and down quark fields respectively. Smeared operators could 
be of value in subsequent studies, particularly for the exploration of the high range. For 
the vector current, j^{x), we use the conserved current, 

jf^i^) = ^u{x){l - -f^)U^{x)u{x + /i) - ^u{x + //)(! + -f^)Ul{x)u{x) (4) 

with fi = 4. In order to extract the matrix element in Eq. from the three-point correlator 
of Eq. (0), one introduces two complete sets of states \n{k)) with the same quantum numbers 
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as in the three-point correlator and gets 



^-(^tf-t)Er„{pf) 



n m 

p-{t-U)E„(p,) 



2Em{pf) 



2E„.{pi) 

This can be simphfied further using 



n{p,)\<P\x)\0) . 



(5) 



(O|0(x)|m(p)) = Z^(p-)e^^-^ 



(6) 



and for a local interpolating field operator, Zm{p) is independent of p. The three-point 
correlator simplifies to 



T^f,^{ti,t,tf,Pi,Pf) =J2J2^r 



{m{Pf)\Jf^WHPi)) 



Z*n- (7) 



Similarly the two-point correlator, which will be needed to get the energies, is given by 



-{t-U)E„{p) 



(8) 



2En{p) 

For periodic boundary conditions on a lattice of Nt time slices, Eq. (jH)) will be modified to 



=E9^e-(t)^"(^-^cosh 



t-U-^]E„{p) 



(9) 



The long time behaviour of the two and three-point correlators will be dominated by con- 
tributions from the lightest pseudoscalar state, i.e. the pion. This asymptotic behaviour is 
given by 

^^~{t~U)E„{pi)-{tf-t)E^{pf) 



r (-1- f f ^ ^ ^ 1 7 



AE^{pi)E^{pf) 



Gtx-k (ti ) t, p) 



|t-ti|>0 



em' 



^-(P^cosh 



t-ti 



(10) 

(11) 



To obtain a reliable result for the pion form factor, we will allow for contributions from 
excited states. For the conserved current, the corresponding transition matrix elements are 
included as follows, 

(p/ - p1) 



{Pf + P 



1 1 11 



(Pf-Pi)' 



-Xpf-Pi)^ 



(12) 
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III. THE ACTION AND ITS PARAMETERS 

Our simulations use the standard Wilson gauge action with P = 6.0. An ensemble of 100 
quenched configurations of size 16^ x 48 was created using a pseudo-heatbath algorithm, 
with 5000 sweeps omitted between saved configurations. The lattice tmQCD fermion action 
is 



SF[^,^P,U]=Y.^P{ 



X] 



1 + 2/€/iZ75r3 + E 7. (V. + V: - V:V,) 



ij{x) (13) 



where the forward and backward lattice derivatives are defined as usual, 



V^ip{x) = U^{x)i;{x + u) -i;{x), (14) 
V*V'(a;) = ^pix) - Ul{x - 0)i,{x - z>), (15) 

and %Ij{x) denotes the doublet of up and down quarks. When ^ = 0, Sp becomes the standard 
Wilson action. Throughout this work, we hold the hopping parameter fixed at its critical 
value, Kc = 0. 156911, [l6| leaving the quark mass directly proportional to /i. Our pion form 
factor studies are carried out with fi = 0.030 and /i = 0.015. We also choose the temporal 
separation of the source and sink to he tf — ti = 15 time steps. The additional cases of 
/i = 0.007, 0.003 and 0.001 are used below as insightful exercises for GMRES-DR. Periodic 
boundary conditions are used in all directions. 

Frezzotti and Rossi have shown that, when the hopping parameter is set to its critical 
value (so-called "maximal twist"), masses and correlation functions with vanishing spatial 

n 

momenta are automatically 0(a) improved in tmQCD. |5| A generic matrix element with 
non-zero spatial momenta can be improved by averaging over momenta of equal magnitude 
but opposite sign as follows. 



'J, k\B\i,pJ + Tjij^B (/, -k\B\i, -pj = 2^b (/, k\B\i,pj Icontinuum + O(a^) (16) 

where ?7j,/,B = ±1 is an overall parity (see Ref. for the precise definition of this parity) 
for the matrix element between the initial state \i,p}, the final state \f,k), and the operator 
B. The renormalization coeffecient relates the continuum and lattice operators. Since 
the energies obtained from a two-point correlator depend only on these energies, like 
masses, are automatically improved without momentum averaging. 
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IV. MATRIX INVERSION 

Some of the standard matrix inverters used in lattice QCD research, such as the stabi- 
hzed biconjugate gradient, fail when applied to tmQCD at maximal twist for sufficiently 
light quarks. |l4l| Fortunately there are other inversion algorithms that succeed for tmQCD 
inversions, such as conjugate gradient, conjugate gradient squared and GMRES.Ql The 
present work made use of the GMRES-DR algorithm^, 

and the remainder of this section 
contains some information about our experience with this inverter. 

The GMRES-DR inverter is built on the standard GMRES (generalized minimal resid- 
ual) matrix inverter, but extends it to incorporate deflation (D) of the smallest eigenvalues 
even after subsequent restarts (R) of the basic GMRES algorithm. Since GMRES-DR is a 
signiflcant improvement over standard GMRES, and since standard GMRES can success- 
fully invert tmQCD matrices, it is interesting to explore the application of GMRES-DR to 
tmQCD. 

GMRES uses a Krylov vector space of some dimensionality (let's call it n) chosen by the 
user and GMRES-DR identifles and retains the /c-dimensional subspace spanned by light 
eigenvectors, where k is also chosen by the user. For the present work, n and k were chosen 
to minimize the wall clock time needed to reach a residual of |r| < 10^^, where r = b — Mx 
for Dirac matrix M and source vector h. This optimization was done at k = and /i = 0.030 
and for our implementation of GMRES-DR(?2,A;) the result was (n, k) = (40, 10). 

For our ensemble of 100 conflgurations, all GMRES-DR(40,10) inversions were successful 
at /i = 0.030, 0.015 and 0.007. The pion form factor was not computed at /i = 0.007 due to 
the onset of flnite volume effects, but the pseudoscalar two-point correlator was computed for 
fi < 0.007 as a means of gaining some experience with GMRES-DR. At yU = 0.003, GMRES- 
DR(40,10) failed to compute one column out of 1200 but increasing the Krylov subspace to 
GMRES-DR(60,10) brought success. At /i = 0.001, GMRES-DR(40,10) failed to compute 
three columns out of 1200 but GMRES-DR(60,10) was again completely successful. Recall 
that our choice of (n, k) = (40, 10) arose from optimization at /i = 0.030; we did not optimize 
separately at these very small values. 

Table H] displays the average number of matrix- vector products that were computed to 
obtain one column of the inverse to a residual of |r| < 10^^ using GMRES-DR(40,10). Since 
this number of matrix-vector products depends on our particular source (i.e. a point source 
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with specific color index and Dirac index) and also on our particular initial value for the 
solution vector, it is more useful to report the change in |r| relative to its initial (i.e. before 
any GMRES-DR iterations) value. In the present case, the initial residual was |ro| = 40.62 
so the data in Table |l] represent the number of matrix- vector products computed to reach 



|r/ro| < 2.5 x 10"®. ' 
other source vectors. 



his is the quantity that can be meaningfully compared to studies with 



17( 1 Figure El shows the pseudoscalar mass squared as a function of fi 



as well as the result of fitting the two largest fi data points to a straight line through the 
origin. Finite volume effects are apparent for /i < 0.007. 



V. RESULTS 



We first analyze the pseudoscalar and vector two-point correlators. Energies are obtained 
by fitting the pseudoscalar correlators to Eq. and the vector correlators to the analogous 
expression. Spatial components of nonzero momenta are averaged over all spatial directions 
to improve statistics; the three spatial components of the vector operator are also averaged. 

Single state fits to the data show a convincing ground state signal for < 4p^ij^, where 
Pmin = 27t/L and L = 16 is the spatial size of our lattice. Multi-state fits were also performed, 
and led to consistent results for the ground state energies. These results can be compared 
to the predictions of the continuum and lattice dispersion relations given by 

(aE,^,y = (aMf + \p\\ (17) 
sinh^(^) =smh^(=f)+i:-^(|). (18) 

respectively. Figures El IH El and El show this comparison where the mass parameters in 
Eqs. ()17|1 and ()18j) were chosen to match the lattice data at p = 0. Table |nl contains the 
numerical values of the pion and p meson masses at those two values of fi for which the pion 
form factor is calculated. 

To extract the pion form factor, we performed a simultaneous fit over the pseudoscalar 
two-point correlator with momentum pi, the pseudoscalar two-point correlator with momen- 
tum pj, and the pseudoscalar (pj)-vector-pseudoscalar(p}-) three-point correlator. The fourth 
component of the conserved vector current was used. To verify the stability of the ground 
state, we've performed both a single state fit over the large time ranges of the two-point 
and three-point correlators where the ground state pion dominates using Eqs. (fTUI) and (fTT|l. 
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and a three state fit over the entire time range (except the source time step) where the 
ground state pion as well as first and second excited states are included. This latter method 
involves 3^ = 9 form factors, Fap-, from Eq. (fT^ . For clarity, here are the explicit forms of 
the correlators used for the three state fit: 



Gl^'^^^iU, t,p,) = E -^e-(^)^"®) cosh \^(t-U-^j EM 



3 



zlzl 



(19) 



G{r'(i.i,P» = E|^e-(^)^'^(^^/)cosh (t-U-^i]EMf) , (20) 



„=1 En{Pf) 



2 



3 3 ^-.{tj-t)Em{pf) ^-{t-U)En{pi) 

{ti, t, tf,pi,pf) — 

n=lm=l ^E^KPf) 2,hn[pi) 

As will be shown, the ground state is quite stable regardless of whether the single state fit 
or three state fit is used. The intermediate case of fitting to a ground state plus one excited 
state leads to similar results, and will not be discussed further. 

For the single state fit, the fitting parameters are the energies the prefactors Z, and 
the form factor -F(Q^) from Eqs. (|iup and (jllj) . For the three state fit, the fitting parameters 
are the energies E and prefactors Z from Eqs. p9|) . ()20|) and ()2H) as well as the form factors 
F{Q'^) and F^p from Eqs. (0) and (fT^ . A standard unconstrained minimization fitting 
procedure was used. Results from a single state fit to the large (Euclidean) time region were 
found to be consistent with results from a three state fit beginning just one time step from 
the source. 

Various choices for pi and pj correspond to comparable momentum transfers, Q^, but the 
cleanest data come when \pi\ and \pf \ are minimal. It is therefore not desirable to restrict 
oneself exclusively to pf = 0. For fx = 0.015 we have computed with pj = (0,0, ±Pmm) as 
well as Pf = 0. For /i = 0.030 we only computed with pf = (0, 0, ±pmin)- Each new value of 
Pf required a new matrix inversion, since the calculation was performed by combining the 
propagator from Xi to x with the sequential propagator from Xi to Xf to x. (See Fig. HI) 
All values of Pi that produce the same q were averaged. Finally, the momentum averaging 
procedure of Eq. ()16|) was employed to remove 0{a) errors. 

In Tables IIIII and IIVI we show the results for the pion form factor obtained from one 
state fits at /i = 0.030 and 0.015 respectively. These same data are displayed graphically 
in Figs. IZIandlHl The physical scale was set to a = 0.10 fm^^. The data show agreement 
with the corresponding vector meson dominance curves. Results from three state fitting are 
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completely consistent with the one state fits, as shown in Tables |V] and IVII and Figs. IHl and 
ITUl Notice the advantage of using a conserved current: the normalization of the form factor 
is F{0) = 1 without any multiplicative renormalization factor Zy- 

In Fig. mj we compare our results at = 0.015 to other quenched calculations of the 
pion form factor at similar quark masses. Two vector meson dominance (VMD) curves are 
also shown — one with the physical p meson mass and the other with the vector meson mass 
taken from our tmQCD simulation at yU = 0.015. The available experimental data are known 
to follow VMD with the physical p meson mass. As evident from Fig. the Wilson results 
have a large systematic lattice discretization error while the tmQCD results are consistent 
with other improved action results and with VMD. Since the quark mass in the Wilson 
simulation of Fig. ^2 is somewhat larger than the quark masses in the other simulations, one 
would expect the Wilson form factor to be slightly larger than the others. The plot clearly 
shows the opposite effect, suggesting that the 0(a) contributions are large. As discussed in 
Ref. yi the apparent smallness of the Wilson form factor is correlated with the smallness of 
the Wilson vector meson mass|19l|. For example, the Wilson form factor would be consistent 
with VMD if the physical scale (i.e. the lattice spacing needed to compute values of 
for the horizontal axis in Fig. [TT|) were obtained from the vector meson mass itself. The 
tmQCD action does not have large 0{a) contributions to the vector meson mass 16| and we 
find correspondingly small 0{a) contributions to the pion form factor. 

As discussed above, the tmQCD results are improved through momentum averaging at 
maximal twist (k = Kc)- Unimproved tmQCD results would be obtained simply by omitting 
the momentum averaging step. Figure IT^ shows the tmQCD pion form factor results aX p = 
0.015 obtained from a one state fit with and without the averaging procedure. Interestingly, 
the data are quite consistent within the statistical uncertainties of our simulation. 
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TABLE I: The average number of matrix-vector products required to compute one column of the 
inverse to a residual of |r| < (2.5 x 10~®)|ro| using GMRES-DR(40,10), where tq denotes the initial 
residual. 



/X average number of products 



0.030 


816 


0.015 


1175 


0.007 


1384 


0.003 


1450 


0.001 


1489 



TABLE II: Ground state pseudoscalar and vector meson masses. 

H 0.015 0.030 
am^ 0.240(3) 0.334(2) 
amp 0.439(27) 0.499(15) 
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TABLE III: Momentum components in units of 27r/L used to extract the pion form factor at 
ji = 0.030 from a single state fit. 0{a) improvement has been invoked by averaging over data with 
n = 1 and n = —1. 





9 [¥] 

^ I J_/ -i 


[GeV2] 

T L J 




(0,0,n) 


(0,0,0) 





1.034 ±0.100 


(0,0,n) 


(0.0. -n) 


0.484 ± 0.025 


0.706 ± 0.056 


(0,0,n) 


(0,±l,-n) 


1.198 ±0.007 


0.454 ± 0.052 


(0,0,ra) 


(±l,±l,-n) 


1.778 ± 0.045 


0.370 ± 0.086 


(0,0,n) 


(0,0,-2n) 


2.402 ± 


0.242 ± 0.034 


(0,0,n) 


(0,±l,-2n) 


2.896 ± 0.1 


0.229 ± 0.044 


(0,0,n) 


(±l,±l,-2n) 


3.281 ± 0.323 


0.327 ± 0.106 


(0,0,n) 


(0,0,-3ra) 


5.242 ± 0.323 


0.161 ±0.065 



TABLE IV: Momentum components in units of 27r/L used to extract the pion form factor at 
/i = 0.015 from a single state fit. 0{a) improvement has been invoked by averaging over data with 
n = 1 and n = —1. 

(0,0,0) (0,0,0) 1.037 ±0.046 

(0,0,0) (0,0,±1) 0.448 ± 0.063 0.669 ± 0.124 
(0,0,n) (0,±l,-n) 1.197 ±0.011 0.401 ±0.110 
(0,0,n) (±l,±l,-n) 1.713 ± 0.059 0.423 ± 0.127 
(0,0,n) (0,0,-2n) 2.401 ± 0.005 0.200 ± 0.066 
(0,0,n) (0,±l,-2n) 2.916 ± 0.056 0.120 ± 0.083 
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TABLE V: Momentum components in units of 2-k/L used to extract the pion form factor at 
/X = 0.030 from a three state fit. 0{a) improvement has been invoked by averaging over data with 
n = 1 and n = —1. 
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Q2 [GeV2] 

T L J 




(0,0,n) 


(0,0,0) 





1.0 ±0.091 




yyj.V). — il ) 


U.^oD in U.Uoo 


n 787 -1- n nciQ 


(0,0,n) 


(0,±l,-n) 


1.197 ±0 


0.473 ± 0.077 


(0,0,ra) 


(±l,±l,-n) 


1.756 ± 0.063 


0.401 ± 0.085 


(0,0,n) 


(0,0,-2n) 


2.394 ± 


0.283 ± 0.041 


(0,0,n) 


(0,±l,-2n) 


2.94 ± 0.042 


0.258 ± 0.055 


(0,0,n) 


(±l,±l,-2n) 


3.294 ± 0.235 


0.37 ±0.124 


(0,0,n) 


(0,0,-3ra) 


5.0 ±0.257 


0.188 ±0.134 



TABLE VI: Momentum components in units of 2it/L used to extract the pion form factor at 
ji = 0.015 from a three state fit. 0{a) improvement has been invoked by averaging over data with 
n = 1 and n = —1. 

(0,0,0) (0,0,0) 0.983 ±0.026 

(0,0,0) (0,0,±1) 0.429 ± 0.055 0.563 ± 0.055 
(0,0,n) (0,±l,-n) 1.197 ±0 0.332 ±0.138 
(0,0,-1) (±l,±l,-n) 1.751 ±0.171 0.44 ±0.218 
(0,0,-1) (0,0,-2n) 2.394 ±0 0.212 ± 0.092 
(0,0,-1) (0,±l,-2n) 2.89 ±0.141 0.115 ±0.104 




anti-quark 



FIG. 1: Three-point correlator for the pion form factor. 
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FIG. 2: Pseudoscalar mass versus the twisted quark mass parameter /x. The solid curve is a fit 
the data at fi = 0.030 and 0.015 to a straight line through the origin. Data points below fi = 0. 
were computed for the exploration of GMRES-DR, not for phenomenological use. 
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FIG. 3: Ground state and first excited state pseudoscalar meson energies from three-state fits to 
tmQCD at /X = 0.030, as compared to the continuum and lattice dispersion relations. 




FIG. 4: Ground state and first excited state pseudoscalar meson energies from three-state fits to 
tmQCD at /Lt = 0.015, as compared to the continuum and lattice dispersion relations. 



16 



2.5 



2- 



1.5- 



0.5- 







T 



T 



T 



• tmQCD, |i=0.030 

— continuum dispersion relation 
□ lattice dispersion relation 



1 excited state 




ground state 



1 



1 







2 3 
p in units of (7U/8) 



FIG. 5: Ground state and first excited state vector meson energies from three-state fits to tmQCD 
at /X = 0.030, as compared to the continuum and lattice dispersion relations. 




FIG. 6: Ground state and first excited state vector meson energies from three-state fits to tmQCD 
at /X = 0.015, as compared to the continuum and lattice dispersion relations. 
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FIG. 7: The form factor at ;u = 0.030 obtained from single state fitting as compared to vector 
meson dominance. 
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FIG. 8: The form factor at = 0.015 obtained from single state fitting as compared to vector 
meson dominance. 
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FIG. 9: The form factor at /x = 0.030 obtained from three state fitting as compared to vector 
meson dominance. 
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FIG. 10: The form factor at fi = 0.015 obtained from three state fitting as compared to vector 
meson dominance. 
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FIG. 11: Comparing the quenched tmQCD pion form factor with other quenched lattice calcu- 

[l(] |. domain wall from Ref. 
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FIG. 12: The pion form factor from tmQCD at yu = 0.015 with and without momentum averaging 



